clc
clear all

load('DP_Data');

x = [0:3500]';
for y=1:length(x)
    FF(y,1) = sum(DP_0c==y-1) / length(DP_0c);
    FF(y,2) = sum(DP_01c==y-1) / length(DP_01c);
    FF(y,3) = sum(DP_1c==y-1) / length(DP_1c);
    FF(y,4) = sum(DP_4c==y-1) / length(DP_4c);
    FF(y,5) = sum(DP_10c==y-1) / length(DP_10c);
    FF(y,6) = sum(DP_40c==y-1) / length(DP_40c);
    FF(y,7) = sum(DP_80c==y-1) / length(DP_80c);
end

mu = [mean(DP_0c) mean(DP_01c) mean(DP_1c) mean(DP_4c) mean(DP_10c) mean(DP_40c) mean(DP_80c)];
wLin = @(xx,alpha) (alpha.*xx);
wBonus = @(xx,B) (B.*(xx>2000));
WW = [wLin(x,0) wLin(x,0.001) wLin(x,0.01) wLin(x,0.04) wLin(x,0.1) wBonus(x,40) wBonus(x,80)];

cnt = DP_Kcombinations;
cnt(1:21,:) = [];

for mm=1:20
m = mm/20;

for i=1:7
    Profit(i) = m*mu(i) - WW(:,i)'*FF(:,i);
end

WeightsRec = zeros(8,length(cnt));
maxProfitRec = zeros(3,length(cnt));

options = optimoptions('simulannealbnd','Display','off');

for counter = 1:length(cnt)
    [m counter]

clear v delta Lstar LfixOpt
KK = cnt(counter,find(cnt(counter,:)>0));
[counter max(Profit(KK))];

K = length(KK);
for k=1:K
    for j=1:K
        v(k,j) = WW(:,KK(k))' * FF(:,KK(j)); % Payment if w(k) is offered and action F(j) is taken
    end
    delta(k,:) = v(k,:) - v(k,k);
end

Linit = exprnd(10,K,K+2);
%Lfix0 = [1 (Profit(KK) ./ mean(Profit(KK))).^50];
%Lfix0 = [1 (10.*Profit(KK) ./ mean(Profit(KK))).^2];
Lfix0 = [0.1 0.1+100.*(Profit(KK)==max(Profit(KK)))];

fun_Lfix = @(Lfix) PrincipalStep1(Linit,Lfix,v,delta,m,mu(KK));
fun_constr = @(LL,Lfix) ConstraintsK(LL,Lfix,v);

[LfixOpt, Objective, exit1] = simulannealbnd(fun_Lfix,Lfix0,zeros(1,K+1),1000.*ones(1,K+1),options);
%[LfixOpt, Objective, exit1] = simulannealbnd(fun_Lfix,Lfix0,-0.01.*ones(1,K+1),[],options);  
%[LfixOpt, Objective, exit1] = fmincon(fun_Lfix,LfixOpt,[],[],[],[],zeros(1,K+1),[])

[obj, Lstar, exit2] = PrincipalStep1(Linit,LfixOpt,v,delta,m,mu(KK));

WeightsRec(1:K+1,counter) = [1; Lstar(1:K,K+2)]./sum(LfixOpt);
maxProfitRec(:,counter) = [-Objective; max(Profit(KK)); -Objective > max(Profit(KK))];
% Line 1: Profit of robustly optimal contract
% Line 2: Profit of the most profitable of the known contracts

end

Mixing = find(maxProfitRec(1,:)>maxProfitRec(2,:));

save(['NewRobustK_' num2str(m) '.mat']);

end